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ABSTRACT 

A new scheme for incorporating radiative cooling in hydrodynamical codes is presented, centered 
around exact integration of the governing semi-discrete cooling equation. Using benchmark calcu- 
lations based on the cooling downstream of a radiative shock, I demonstrate that the new scheme 
outperforms traditional explicit and implicit approaches in terms of accuracy, while remaining com- 
petitive in terms of execution speed. 

Subject headings: hydrodynamics — methods: numerical — radiation mechanisms: thermal — shock 



1. INTRODUCTION 

Cooling by optically thin radiative emission plays 
an important role in many differing types of astro- 
physical flow. This is especially true of radiative 
shocks, which arise in a wide variety of co ntexts (e.g., 
iStrickland fc Blondinl fl995t iMignond [20051 . and refer- 
ences therein). In this paper, I present a new scheme for 
incorporating radiative cooling in hydrodynamical codes, 
centered around exact integration of the governing cool- 
ing equation. 

To lay the necessary groundwork for the subsequent 
discussion, Section [2] presents a derivation of the semi- 
and fully discrete forms for the cooling equation. Sec- 
tion [3] then reviews various schemes used to solve this 
equation, culminating with the introduction of the new 
exact integration (EI) scheme. These schemes are bench- 
marked in Section |4] to explore the relative trade-offs be- 
tween accuracy and execution speed, and I conclude with 
brief remarks in Section [5j 

2. THE COOLING EQUATION 

The hydrodynamical equation of energy conservation 
for an ideal gas can be written as 



-( 7 -l)n e n H A(T) 



(1) 



dP jP dp 
~dt ~p~dt 

Here, P is the pressure, p the density, T the temperature, 
7 the ratio of specific heats, n c and nu the electron and 
hydrogen number densities, respectively, and d/dt de- 
notes the Lagrangian (total) time derivative. The func- 
tion A(T) represents the electron cooling efficiency, and 
is typically o btained in tabular form from detailed model- 
ing (see. e.g..lRavmond et al.lll976HSutherland fc Dopital 
Il99a lOnat fc Sternberg! 12007ft The dependence of A 
on temperature alone is ultimately what makes the EI 
scheme possible, but is also somewhat of an idealization 
of the underlying physics. More-sophisticated treatments 
incorporate additional explicit dependencies on ioniza- 
tion balance, by tracki ng a time- varying network of ionic 
abun dances (see, e.g.. iRaga et al.l 12000: Mignone et alj 
12007ft . It is not yet clear how the EI scheme might be 
extended to these treatments. 

In Eulerian-based, finite-difference hydrodynamic 
codes, it is common to implement energy conservation (jTJ) 
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using an operator splitting approach. The rate of pres- 
sure change is divided into a component associated with 
adiabatic expansion/contraction 



dP _ jP dp 
dt ad p di' 

and a component associated with radiative cooling, 



dP 
~dt 



-(7 - l)n e n H A(T). 



(2) 



(3) 



cool 



The pressure change due to the adiabatic component |(2]) 
is typically applied in the advection stage of the code, 
during which the density and velocity are also updated 
in accordance with the mass and momentum conserva- 
tion equations 1 . The pressure change due to the cool- 
ing component ([3|) is then applied in a subsequ ent stage, 
durin g which the density is held constant (e.g.. IMignond 
12005ft . The isochoric nature of this latter stage does not 
preclude si mulations of isobaric syst ems such as cooling 
flows (e.g.. iPeterson fc Fabian! [2006): in these cases, the 
components (|2l3j) are equal and opposite, resulting in no 
net pressure variations. 
The ideal gas law 



P = nkT, 



(4) 



is used to recast the cooling equation (|3|) in terms of 
temperature; here, n is the total number density of par- 
ticles, and k is Boltzmann's constant. The mean molec- 
ular weight p, = p/n is a ssumed to remain constant; as 
(|Gnat fc Sternberg! 12007! ) argue, this is a reasonable ap- 
proximation for temperatures > 10 4 K (although it can 
break down in circumstances where depa rtures from ion- 
izatio n equilibrium are significant; see Tcsilcanu et alJ 
2008;) . The cooling equation then becomes 



dT 
dt 



(7 ~ 



A(T), 



(5) 



where p c = p/n c and pa = pjn^ are the effective molec- 
ular weights per electron and per hydrogen atom/ion. In 

1 By maintaining an adiabatic advection stage, it remain s pos - 
sible to use numerical schemes derived from Godunov's (|1959D 
characteristic-based approa ch, such as the popular Piecewise 
Parabolic Method (PPM) of lColella fc Woodward! JT984T) 
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the regime of full ionization, the molecular weights ap- 
pearing in this expression are given by 



M = 



Me = 



2Jf + 3(l-X-Z)/4 + Z/2 
2u 



1 + X' 



(6) 

(7) 
(8) 



here, X and Z are the usual hydrogen and metal mass 
fractions, and u is the atomic mass unit. 

Implementing eqn. (J5j) in a hydrodynamical code re- 
quires discretization in space and perhaps also time. For 
simplicity, I focus here on a zeroth-order finite-volume 
spatial discretization 2 , which leads to the semi-discrete 
cooling equation 



d7^ 
dt 



(7 - 1 w 



A(3i) 



(9) 



Here, pi is the volume-averaged density in the numerical 
zone with integer index i, while T{ is the corresponding 
zone temperature, _calculated from pi and the volume- 
averaged pressure Pi using the ideal-gas law ((4|. 

The following sections discuss various approaches to 
solving this equation across a discrete time step At. The 
explicit fi )3.ip and implicit f ^3.2|) schemes are all based 
on a fully discrete cooling equation, derived from the 
semi-discrete form by replacing the temperature rate- 
of-change with a finite difference: 



T 



n+l 



rpn 



(7 - 1 w 



At 



ACT,). 



(10) 



The superscripts n and n + l indicate values at consec- 
utive times t" and t n+1 = t n + At. Because the finite 
difference is centered, this equation is second-order ac- 
curate in time. To evaluate the cooling efficiency A(Ti) 
on the right-hand side, either the initial or the updated 
temperature may be used; the choice differentiates ex- 
plicit schemes from implicit schemes. 

3. SOLVING THE COOLING EQUATION 

3.1. Explicit Schemes 

In an explicit scheme, the cooling efficiency in eqn. (I10| 
is evaluated using the initial temperature T": 

nn+l rpn 



rprt 



(7 - 1 w 



A(37). 



(11) 



At ku e (j, H 

This is now first-order accurate in time, because the 
right-hand side is not centered in the interval (t™, t n+1 ). 
Solving for the updated temperature, 



where 



rpn- 



^cool — 



rpn 



1 



At 

^cool 



( 7 - l)p 4 /iA(T«) 



(12) 



(13) 



is the single-point cooling time. This solution, together 
with the ideal gas law ((4]), allows the pressure in each 
zone to be updated across the time step At. 



2 Sec Strickland & Blondin (1995) for a demonstration of how a 
higher-order discretization might be constructed. 



The behavior of the explicit scheme ([12]) is investigated 
by calculating the updated temperature f" + as a func- 
tion of time step, for three different choices of the initial 
temperature: T? = (10 6 K, 10 7 K, 10 8 K). The cooling ef- 
ficiency is obtained from a piecewise power-law fit to the 
col lisional ionization equilibr ium (CIE) values tabulated 
by iGnat fc Sternberg! {2007) . Because the tabulation is 
truncated at 10 4 K, this temperature is imposed as floor 
on y™ +1 (in an actual simulation, this floor temperature 
might correspond to the reheating effects of a nearby 
star). A monatomic gas (7 = 5/3) and solar abundances 
(X = 0.7, Z = 0.02) are assumed here and through- 
out. The upper panels of Fig. [1] plot the results from 
these calculations. By way of comparison, the panels 
also show the exact solutions to the semi-discrete cool- 
ing equation ([9]); Section [3~3l discusses how these solu- 
tions are obtained. 

For At approaching t CO oi, the explicit scheme leads 
to updated temperatures that depart quite significantly 
from the exact values. Put simply, this is because the 
cooling efficiency is fixed at its initial value A(T"), rather 
than being allowed to evolve in response to the cooling 
process. This difficulty can be avoided by dividing the 
time step At (typically set during the advection stage by 
the Courant-Friedrichs-Lewy criterio n; see $H) into a se- 
quen ce of smaller sub-steps (see, e.g.. lPlewa fc Rozvczkal 
1992), and applying eqn. (TT2"]) multiple times. Alter- 
natively, a higher-order temporal discretization of the 
cooling equation is possible; for instance, a second-order 
Runge-Kutta method has 



T. 



n+1/2 



rpn 



1 At 



2t 



cool 



1 



A(r" +1/2 ) At 

A(T") t cool 



(14) 
(15) 



The lower panels of Fig. Q] plot the updated tempera- 
tures calculated using this second-order scheme. There 
is a clear improvement over the first-order approach, and 
further improvemen ts can be gained by go ing to even- 
higher orders (e.g.. iSutherland et al.l [2003h . However, 
with each order added an additional evaluation of the 
cooling efficiency is required; hence, the computational 
costs necessarily escalate. 



To overcome 
when At > t coo i , 
Strickland & Blondin 



3.2. Implicit Schemes 
the drawbacks of explicit schemes 



a nu mbe r of autho rs (e.g. , 
fl995t iStone et alJ Il997t 
Pittard et al.ll2004f ) instead opt for an implicit scheme. 
The cooling efficiency in eqn. (jTOJ) is then evaluated 
using the updated temperature T" + : 

T z n+1 (7 ~ 



A(7f +1 ). 



(16) 



At fc/Z /iH 

The solution can be written in a standard form similar 
to the explicit case (cf. eqn. [T2"]) . 

A(If +1 ) At 



^n+l 



rpn 



1 acH») t cool J' (17) 

but the appearance of r™ +1 on the right-hand side means 
that this equation must now be solved numerically, typ- 
ically using a root-finding algorithm. 



An Exact Integration Scheme for Radiative Cooling in Hydrodynamical Simulations 



3 




.0 0.1 1.0 0.1 1.0 

(Co,) A« (i cool ) At (Z cool ) 

Fig. 1. — The updated temperature TV" plotted as a function of time step At (in units of the cooling time t coo i), for three differing 
choices (left-to-right) of initial temperature T™. The circles in the top (bottom) panels indicate values calculated using the first-order 
(second-order) explicit scheme. The solid lines show the corresponding exact solutions. 



T t n = 10 6 K 



f, n = 10 7 K 




Ai (C=i) 



A« (Col) 



At (Col) 



Fig. 2. — As in Fig.[TJ except that the circles in the top (bottom) panels now indicate values calculated using the secant (Brent) first-order 
implicit scheme. 



Fig. [5] investigates the behavior of this implicit scheme, 
plotting TJ 1 as a function of At for the same parame- 
ters as in Fig. [TJ The upper panels use a secant algorithm 
to solve eqn. (TTT]) , with a fallback to bisection when the 
most recent iterate for J" ™ +1 falls outside the bounds of 
the A(T) tabulati on. Conversely, t he lower panels use 
Brent's algorithm (Pre ss et al.| [T992). In both cases, so- 
lutions are iterated until the fractional change in T" +1 
drops below 10~ 4 . 

The figure reveals problems with the implicit schemes. 
For instance, in the T™ = 10 6 K case, T™ +1 i s signifi- 
cantly underestimated for 0.3i coo i ^ At < 0.7i coo i, and 



overestimated for At > 0.7t coo \. Moreover, rather than 
varying smoothly as At is increased (as one might hope 
for a stable scheme), T" +1 exhibits abrupt jumps. 

To explore the origin of these jumps, I introduce the 
twin discriminants 

=.«+! 



and 



D a =T?- 

„A(T/ l+1 ) At 



Dh = T: 



A(I7*) icooi' 



(18) 
(19) 



such that the implicit equation (fT7|) corresponds to the 
condition D a = Df,. Fig. [3] plots the discriminants to- 
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At = 0.1 t 



At = 0.3 / 



At = 2.0 £„ 
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Fig. 3. — The discriminants D a (thick) and D b (thin) plotted as a function of updated temperature T™ +1 , for three differing choices 
(left-to-right) of the time step At. Intersections between the two curves are highlighted by circles. 
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Fig. 4. — As in Fig. [T] except that the circles in the top (bottom) panels now indicate values calculated using the secant (Brent) 
second-order implicit scheme. 



gether as a function of T? +1 , for T™ = 10 6 K and three 
choices of time step. The middle, At = 0.3i coo i panel 
shows that the curves intersect multiple times, corre- 
sponding to multiple, distinct solutions (in this case, five) 
to the implicit equation. Abrupt switching between these 
solutions, due both to the convergence behavior of the 
particular root-finding algorithm, and to the appearance 
or disappearance of solutions as At is varied, is respon- 
sible for the jumps seen in Fig. O 

To underscore further that solution jumping is an in- 
trinsic property of implicit schemes, Fig. Q] illustrates so- 
lutions of the cooling equation using the Crank-Nicholson 
method, 



rpn+l rpr, 

At 



kflefJ-H 2 



(20) 



which is now second-order implicit. In the standard form, 
this becomes 



rj-vn 



l 



A(T™ +i ) + A(77) At 



2 A (77' 



^cool 



(21) 



While the data plotted in the figure differ from the first- 
order cases shown in Fig. [2l they still exhibit abrupt 



jumps arising from the existence of multiple solutions 
— although the range of At values over which jumping 
occurs is somewhat reduced. 

In spite of these various issues, implicit schemes 
have proven popular in the literature. This stems 
in part from their reputati o n for stability; for in- 
stance, IStrickland &; Blondinl |l995) remark that their 
implicit cooling scheme 'i s unconditionally sta ble' (words 
subsequently echo ed bv iPittard et al.1 12004): likewise, 
I Stone et"aTl (|1997f l state that their scheme l is stable even 
when the cooling time is much less than the dynami- 
cal time\ However, this confidence appears misplaced. 
While implicit schemes are stable when used to solve lin- 
ear equations (as can be demonstra ted through a vo n 
Neumann stability analysis; see, e.g. JPress et al.l Fl992). 
this property does not necessarily extend to non-linear 
systems such as the semi-discrete cooling equation dH). 
Moreover, even if stability can be established for a given 
scheme, there are no corresponding guarantees of accu- 
racy or convergence — and it is in these latter capacities 
that the implicit schemes reviewed here fall short. The 
jumping between solutions is particularly problematic, 
because it can cause two neighboring zones with very 
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similar initial states to cool to quite different temper- 
atures. This will establish a strong pressure differential 
between the zones, in turn generating spurious fluid flows 
and/or waves. 

3.3. The Exact Integration Scheme 

The new cooling scheme introduced here avoids the 
various difficulties outlined above, by going back to the 
semi-discrete cooling equation (|9]) and solving it exactly. 
The equation is first rearranged as 

dT_i = (7 ~ l )Pi^ A t 
A(Ti) A:/i G ^ H 

and then integrated across a time step: 

r T ? +1 _w_ = (7 - 



At. 



(22) 



(23) 



The dimensionless 'temporal evolution function' (TEF) 
is then introduced as 



Y(T) 



A(T ref ) 



dT' 



(24) 



r ref j T A(r') 

where T re { is an arbitrary reference temperature; the 
TEF represents a normalized measure of the total time 
taken to cool from T re f to T. With this definition, the 
integrated cooling equation ([25)1 becomes 



[Y(TP ) - Y(T?^)] = 



At 



(25) 



T? A(T ref ) 

where t coo \ has the same definition as before, and the 
consequent solution is 



rpn+l _ y- 



Y(fP 



f? A(T ref ) At 
T rcl A(T") tcooi 



(26) 



This result is exact, but requires construction of the TEF 
and its inverse from A(T). The Appendix presents an- 
alytic expressions for Y(T) and its inverse in the com- 
mon cases where the cooling efficiency is represented by a 
power law ( flA.l)) and a piecewise power law ( §A.2[) . The 
cooling efficiencies used in Figs.Q]and[2]fall into the latter 
category, and the exact solutions plotted in these figures 
are calculated using the EI scheme desc ribed here . Like - 
wise, the cooling efficiency assumed by iMignond (|2005[ ) 
falls into the former category, and in fact his analytic 
cooling scheme (which foreshadows the present paper) 
can be derived from this EI formalism. 

4. BENCHMARKS 

The preceding sections (and in particular, Figs. [T] 
and [2]) demonstrate that both explicit and implicit 
schemes for solving the cooling equation can become in- 
accurate as the time step approaches the cooling time 
£ C ooi; in contrast, the EI scheme gives the exact solution 
for any value of At. However, an important caveat here 
is that the time step is itself constrained by numerical 
considerations in the advection stage. Efficiency dictates 
that At be chosen as large as possible (subject to ac- 
curacy requirements), but for stability reasons it cannot 
exceed the limit established by the Couran t-Friedrichs- 
Lewy (CFL) criterion (see, e.g.. lLanevlll998[ ). 

To explore how the differing cooling schemes perform 
with a CFL-based time step, I consider the problem of a 



steady, 3 1-dimensional radiative shock characterized by 
an upstream density p ln , Mach number A^m and temper- 
ature Tj n . For various combinations of these parameters 
(to be discussed below), each scheme is benchmarked by 
repeating the following steps: 

1. The run of density and pressure throughout the 
post-shock cooling re gion are calculated using the 
approach described bv lStrickland fc Blondinl (|1995l 
their §4.1). This region is bounded on the upstream 
side by the shock itself, and on the downstream side 
by the condition T = T- m (i.e., the gas has cooled 
back down to its initial temperature). 

2. The cooling region is discretized into N equal-sized 
zones; for each zone, the volume-averaged density 
pi and pressure Pi are evaluated, and the corre- 
sponding temperature Ti is calculated using the 
ideal-gas law ([4]). 

3. The CFL time step is calculated as At = Ax/c max , 
where Ax is the spatial extent of the zones, and 
c max = max(-\/7 I\j pi) is the maximum value of 
the adiabatic sound speed over all zones composing 
the cooling region. 

4. For each zone, the updated temperature T™ +1 is 
evaluated using one of the cooling schemes dis- 
cussed in the preceding sections. This step is re- 
peated 5 times, and the average CPU execution 
time r (per zone, per repeat) is recorded. 

5. The updated temperatures are compared with the 
exact values T"^ that result from using the EI 



scheme; the maximum relative error 

£ = max(|Tf +1 -^+ 1 |/^ 
is recorded. 



EI > 



(27) 



To cover a representative region of parameter space, I 
consider three values M. m = 3, 10, 100 of the Mach num- 
ber (corresponding to mild, moderate and strong shocks), 
and three values N = 1, 10, 100 of the zone count (corre- 
sponding to poor, moderate and good resolution of the 
shocks). An upstream density p m = 10~ 15 gcm -3 is as- 
sumed throughout, but results do not depend at all on 
this value. With these choices of parameters, and for 
each of the five cooling schemes considered previously, 
Table Q] shows the error e and execution time r obtained 
by following the steps above. All calculations were un- 
dertaken on a single core of an Intel E5345 quad-core 
CPU running at 2.33 GHz. 

The table reveals a general trend that the error de- 
creases as the zone count increases. When the shock is 
resolved by only a single zone, the error tends to be large 
(with the obvious exception of the EI scheme, for which 
e = always). For N = 10, e is below 10% in all but 
one case; and by N — 100, it is below 1% in all cases. To 
explain this trend, the definition of the CFL time step is 
used to write 

r 1 - (28 > 

^cool Cmax^cool 

3 In reality, radiative shocks are often time-varia ble due to the 
cooling instability discovered bylLangcr ct al. (19811); however, this 
variability is ignored here since the principal criterion is a well- 
defined test system, even if it is somewhat idealized. 
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TABLE 1 

Benchmark results 





N 


= 1 


N = 


10 


N = 


100 


Cooling Scheme 


e{%) 


r(ns) 


e(%) 


r(ns) 


e{%) 


r(ns) 




M in = 3 


1 st -order explicit 


4.0 


95 


14.1 


111 


0.2 


100 


2 nd -order explicit 


4.0 


213 


5.9 


201 


0.0 


181 


l at -order implicit (secant) 


26.1 


844 


4.2 


484 


0.2 


362 


l st -order implicit (Brent) 


26.1 


903 


4.2 


787 


0.2 


562 


2 nd -order implicit (secant) 


22.0 


931 


2.1 


837 


0.0 


612 


2 nd -order implicit (Brent) 


22.0 


927 


2.1 


837 


0.0 


613 


Exact 


0.0 


213 


0.0 


192 


0.0 


173 



M in = 10 



1 st -order explicit 


30.8 


94 


7.6 


112 


0.3 


101 


2 nd -order explicit 


3.0 


215 


0,1 


204 


0.0 


183 


l st -order implicit (secant) 


24.0 


906 


6.5 


548 


0.2 


383 


l st -order implicit (Brent) 


24.0 


972 


6.5 


778 


0.2 


572 


2 nd -order implicit (secant) 


12.2 


1230 


1.3 


875 


0.0 


640 


2 nd -order implicit (Brent) 


12.2 


1248 


1.3 


874 


0.0 


639 


Exact 


0.0 


215 


0.0 


197 


0.0 


173 



M in = 100 



1 st -order explicit 


38.0 


92 


1.8 


107 


0.1 


98 


2 nd -order explicit 


12.1 


207 


0.1 


204 


0.0 


181 


l st -order implicit (secant) 


99.9 


266 


1.2 


523 


0.1 


387 


l st -order implicit (Brent) 


99.9 


258 


1.2 


725 


0.1 


544 


2 nd -order implicit (secant) 


99.9 


315 


0.3 


796 


0.0 


598 


2 nd -order implicit (Brent) 


99.9 


317 


0.3 


799 


0.0 


600 


Exact 


0.0 


214 


0.0 


195 


0.0 


169 



Because the flow downstream of the radiative shock is 
sub-sonic, it follows that 

^ < (29) 

t-cool ^t-cool 

where v < c max is the typical flow velocity in the cooling 
region. Recognizing that vt coo \ approximates the spatial 
extent of this region, the corollary is that 
At 1 



, ~- v (30) 

tcool JV 

Thus, the limit N 3> 1 implies that At <C i C ooi, which fa- 
vors accurate cooling irrespective of the choice of scheme. 
Turning this statement around, all of the cooling schemes 
apart from the EI scheme tend to be inaccurate when the 
cooling region is poorly resolved. Of course, this applies 
only when At is tied solely to the CFL time step. It is of- 
ten desirable to place further constraints on At, over and 
beyond that given by the CFL criterion. For instance, to 
improve the coupling between thermal and hydrodynam- 
ical evolution, At can be limited so that the anticipated 
temperature/pressure change during cooling does not for 
any zone exceed a specified fraction of its initial value. 
The price paid in this approach is the greater number of 
steps that must be taken to cross a given time interval. 

Looking now at the relative performance of the differ- 
ing schemes, the r data in Table [T] reveal that the first- 
order explicit scheme is the fastest, with an average ex- 
ecution time of ~ 100 ns per zone. The second-order ex- 
plicit scheme and the EI scheme are both only about fac- 
tor of two slower than this, with the latter slightly beat- 
ing the former in all but two tests. The implicit schemes 
are in every case the slowest, ranging from around 2.5 up 
to 12 times slower than the first-order explicit scheme. 



5. CONCLUDING REMARKS 

A criticism that might be leveled at the exact inte- 
gration scheme is that it requires the reciprocal of the 
cooling efficiency, 1/A(T), be analytically integrable. In 
practice, this is rarely an issue; the most common repre- 
sentations of A(T) are piecewise power-law or piecewise 
polynomial fits to detailed models, both of which meet 
this restriction. In any case, it is always possible to re-fit 
arbitrary cooling efficiency data with a conforming rep- 
resentation. 

The principal strengths of the EI scheme are twofold. 
On the one hand, it produces exact solutions to the semi- 
discrete cooling equation Q, irrespective of whether 
the time step is small or large compared to the cool- 
ing time t coo \. On the other, it remains very competi- 
tive in terms of execution speed, being only two times 
slower than the (fastest, yet often inaccurate) first-order 
explicit scheme. While more-sophisticate d cooling treat- 
ment s that track ionic abundances (e.g., Mignone et al. 
2007) will remain the state-of-the-art in terms of physical 
fidelity, the strengths of the EI scheme naturally recom- 
mend it as the cooling scheme of choice in any hydrody- 
namical code where a simple, fast and robust treatment 
of optically thin radiative losses is desired. 
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APPENDIX 

TEMPORAL EVOLUTION FUNCTIONS 

Power Law 

I first consider the simple case of a power-law cooling function, 

T 



A(T) = A ref — 

, J ref 



(Al) 



where T TC { is the reference temperature introduced in N3.31 and A re f and a are constant coefficients. Substituting this 
into eqn. J24|) leads to a TEF 



Y(T) = < 



a=l. 



The corresponding inverse TEF is given by 

Y~\Y) = 



T ref [1 - (1 - a)Yf 
T rcf exp(-Y) 



a = 1. 



(A2) 



(A3) 



Piecewise Power Law 

Mo re physically realistic cooling functions are often represented by piecewise power-law fits to detailed models 
(e.g.. iWalder fc Folinil [l99l iKimoto fc Chernofl Il997t iCaunt fc Korpil l200lt iTownsend et al.l [2007h . I assume a fit 
parametrization of the form 

A(T) = A k (^j T k <T<T k+1 , (A4) 

for a set of N — 1 temperature intervals (Tf., T k +i) (k = 1, 2, . . . , N — 1) and coefficient pairs {A k , a k }- Substituting 
this into eqn. with a reference temperature chosen as T re f = T/v, leads to the piecewise TEF 



1 Ajv T fc 



Y(T)=Y k + { \-<*!l a <° T » 



"fc = 1 



(A5) 



where Aat = Aat_i(T/v/T/v_i) q;jv - 1 . The coefficients Yfe = Y(Tfc) are constants of integration; the requirement that 
F(T) be continuous dictates that 



Yk — Y k+1 — < 



1 Aw T k 
1— otk A fc T N 



1 - -ii. 



A fc Tn 111 VT fc+1 y 



(A6) 



This recurrence can be started by noting that Yjy = V(T rc f) = (cf. eqn. 
(Y k , Y k+1 ) interval by 



afc = 1. 

. The inverse TEF is given in each 



Y~\Y) 




l-(l-c*k)t?ff(Y-Y k ) 



!/(!-£»*) 



«fc 7^ 1 

01/8=1 



r fe < y < y fc+ i. 



(A7) 
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